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ABSTRACT 

Binary population synthesis calculations and associated predictions, especially event rates, are known 
to depend on a significant number of input model parameters with different degrees of sensitivity. At 
the same time, for systems with relatively low formation rates, such simulations are heavily compu- 
tationally demanding and therefore the needed explorations of the high-dimensional parameter space 
require major - often prohibitive - computational resources. In the present study, to better understand 
several key event rates involving binary evolution and binaries with two compact objects in Milky Way- 
like galaxies and to provide ways of reducing the computational costs of complete parameter space 
explorations: (i) we perform a methodical parameter study of the StarTrack population synthesis 
code ; and (ii) we develop a formalism and methodology for the derivation of multi- dimensional jits 
for event rates. We significantly generalize our earlier study, and we focus on ways of thoroughly 
assessing the accuracy of the fits. We anticipate that the efficient tools developed here can be applied 
in lieu of large-scale population calculations and will facilitate the exploration of the dependence of 
rate predictions on a wide range binary evolution parameters. Such explorations can then allow the 
derivation of constraints on these parameters, given empirical rate constraints and accounting for 
fitting errors. Here we describe in detail the principles and practice behind constructing these fits, 
estimating their accuracy, and comparing them with observations in a manner that accounts for their 
errors. 

Subject headings: binaries:close — stars:evolution — stars:neutron - black hole physics 



1. INTRODUCTION 

Models for binary stellar evolution and population syn- 
theses are necessary to provide quantitative theoretical 
predictions for the relative likelihood of assorted events 
involving the evolution of binary stars. The resulting 
predictions are particularly critical when no empirical 
estimates exist for topics of immediate astrophysical in- 
terest, such as mergers of double-compact object (DCO) 
through the emission of gravitational waves. The most 
practical and widely applied binary population synthesis 
codes available in the community - such as t he Star- 
Track code described in Belczynski ct al. (2002) and sig- 
nificantly updated in Belczynski et al., ( .2006a ) : the BSE 
code, descri bed in iHurlev et al.l (l2002f); the S cBa code 
described in lPortegies Zwart fc VerbuniJ (Il996f) : and the 
StarFaster code described in lFrver et aU ()1998[ ) - rely on 
a large set of fairly simple parameterized rules to char- 
acterize many complex and often ill-understood physical 
processes. Unfortunately, current binary population syn- 
thesis codes are greatly and often forbiddingly computa- 
tionally demanding (depending on their level of sophisti- 
cation): even with substantial simplifications, exploring 
the entire parameter space relevant to the simulations is 
beyond present-day computational capability. 

However, observational information can provide us 
with constraints that help us improve our understand- 
ing of massive binary evolution. For example, pulsar 
searches continue to discover and refine observations of 
isolated pulsar s and new binary pulsar systems; e.g., see 
iLorimerl ()2005f ) for a review. Specifically, the samples of 
binary pulsars with neutron star and relatively massive 
white dwarf companions have been used for a statistical 
derivation of empir ical rate es t imate s for their formatio n 
(most recently see iKim et aH (|2003l ). iKim et al.l (|2004j ). 



IKalogera et al.l ()2005| ) and references therein). Addi- 
tionally, many ground-based gravitational wave detectors 
now operating at or near design sensitivity (i.e., LIGO, 
GEO, TAMA) are designed to detect the late stages 
of double compact object (DCO) inspiral and merger. 
Based only on early-stage data, these instruments have 
already provided conserv ative upper limits to certain 
DCO merger rates (see.e.g. [Abbott et al.ll2005bl lal). With 
LIGO now very close to design sensitivity, a year of 
LIGO data could definitively exclude (and even possi- 
bly confirm) the most optimistic t heoretical predictions 
for BH -BH merger rates (see, e.g. lO'Shaughnessv et al.l 
l2005al for a range of BH-BH merger rates arising from 
binary evolution in Milky Way- like galaxies). Thus, 
gravitational-wave based upper limits (and, eventually, 
detections) will shortly provide constraints on theoreti- 
cal models of DCO formation. 

Faced with the availability of empirical rate con- 
straints, and yet at first unable to quantitatively impose 
them on population synthesis predictions because of the 
extremely limite d exploration of the multi-di menional 
parameter space, lO'Shaughnessv et al.l (|2005af l realized 
that any single unambiguous population synthesis pre- 
diction could be sampled loosely and then fit over the 
most sensitive population synthesis parameters. In the 
same study, we also presented a technique to acceler- 
ate the synthesis code used [StarTrack) to study a sin- 
gle target subpopulation, which we called "partitioning" : 
we used experience gained from prior simulations to re- 
ject binary parameters h ighly unlikely to produce the 
current event of interest. lO'Shaughnessv et al.l (|2005bD 
first applied these early fits to allow a direct compar- 
ison between S'torT'racfc-produced population synthesis 
predictions and the observed formation rate for NS-NS 
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binaries. Though only a smaU fraction (2%) of Star- 
Track models appeared consistent with the constraints, 
conceptual challen ges with seven-dimension al visualiza- 
tion prevented O'S haughnessv et al.l ()2005aD from clearly 
describing the constraint-satisfying reg ion. A forthcom- 
ing paper. lO^'Shaughnessv et all (|2007b( l will significantly 
extend that earlier preliminary analysis, adding signifi- 
cantly more observational constraints as well as a clearer 
investigation of the constraint-satisfying models. 

In th is pap er, we extend and generalize the analysis 
of lO'Shaughnessv et al.l (j2005al), and we present a thor- 
ough discussion of our much-updated and vastly larger 
population synthesis archive (Sj2]) and particularly of the 
fitting methods we employ to extract predictions and as- 
sess their uncertainties and quality (^. Because sev- 
eral imporant and yet not immediately obvious pitfalls 
must be identified avoided when constructing, testing, 
and applying fits, in this paper we provide a thorough 
and pedagogical discussion of our fitting method and phi- 
losophy. For example, we explain how systematic fittin; 
errors that were ignored in lO'Shaughnessv et al.l (|2005 
can limit our ability to constrain t he merger rates in the 
Milky Way. In a companion paper ([O'Shaughnessv et al.l 
l2007bf) we explore the astrophysical applications of these 
fits, emphasizing their use in deriving empirical con- 
straints on population syntheses and rate predictions. 

Though we develop the fitting formulation specifically 
for the results of the StarTrack code, the methods de- 
scribed in this paper can be applied to any family of pop- 
ulation synthesis simulations. For this reason, we survey 
the physics underlyin g the StarTrack model o nly in our 
companion paper lO'Shaughnessv et all(|2007bD . in which 
these fits are used to cover the parameter space and apply 
empirical rate constraints from supernovae and binary 
pulsars. 

2. POPULATION SYNTHESIS ARCHIVES 

Population synthesis simulations can be extremely 
computationally demanding: even though StarTrack can 
fully evolve roughly 10'^ binaries of interest^ per CPU- 
hour with modern-day processors, because some double 
compact objects form very infrequently (e.g., black holes, 
which occur roughly once every ^ I0~'* binaries evolved), 
a representative sample of stellar systems often contains 
104-5_io6-5 binaries and requires hundreds of CPU hours 
to complete. Additionally, since population synthesis 
rate predictions depend delicately on model parameters, 
the computation time needed to build up a sufficiently 
representative collection of stellar systems - one where 
some event of interest occurs many times - varies consid- 
erably depending on astrophysical assumptions. Given 
the prohibitive computational demands of a brute-force 
approach, we took advantage of several simplifications 
originally developed in O'Shaughncssy et al. (2005a) to 
assemble our archive of roughly 3000 population synthe- 
sis simulations, upon which our fits our predicated. In 
this section we briefly describe how those archives were 
generated and how we identify and extract event rates 
for several processes of interest. 

2.1. StarTrack population synthesis code 
^ Specifically mi > 4Mq ; see the discussion below. 



We estimate formation and merger rates for several 
classes of double compa ct objects using t he St arTrack 
code first developed by iBelcz vnski et al.l (|2002f ) [here- 
after BKB] and recently significantly upda t ed and tested 
as described in detail in iBelczynski et al.l ([2006a) . Like 
other population synthesis codes, StarTrack evolves some 
number N binaries from their birth (drawn randomly 
from specified birth distributions) to the present, track- 
ing the stellar and binary parameters. Because binaries 
without any high-mass components cannot produce black 
holes or neutron stars, we configured StarTrack to simu- 
late only those binaries whose heaviest initial component 
mass TOi was greater than AMq. In effect, the small num- 
ber N of binaries simulated mimic the results of a much 
larger simulation of size N^ff in which mi can take on 
any value from the hydrogren burning limit (O.OSMq) to 
the maximum initial stellar mass we allow (ISOM©): 
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based on a broken Kroupa initial mass function (j){m) cx 
if m e [0.08,0.5]A/o, oc m'^-^ if m g [0.5,1]Mq, 
and (X m"^'^ if to > IMq. In turn, a simulation of 
Neff stellar systems can be scaled up to represent the 
Ng stellar systems in the Milky Way. The number of 
stellar systems iV^ in the Milky Way (and thus the nor- 
malization of population synthesis simulations) can be 
estimated in many ways, depending on the observational 
inputs used; for the purposes of this paper, we estimate 
Ng from the average number of binary systems formed 
through steady star formation over T = 10 Gyr of star 
form ation at M ~ S.5Mq yr"! ([Ranalll991l : lLacev fc FalU 
[HH): 
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^ (TO1+TO2) (toi) (1 + /b (m2/mi)) 



(2) 



:7.6 X 10^ 



^lOGyrj 1 + /fc (TO2) / (toi) 



where fb is the binary fraction and (TO2) / (toi) is the 
ratio of average masses of the companion and primary, 
respectively. The proportionality constant s = Ng/Neff 
needed to scale the results of our simulations up to the 
universe is therefore 



S^Ng/N^ff 

= 5.9 X lG*^iV~i(T/10Gyr) 



(3) 



l + /fc[(m2)/(TOi)] 



Note that parameters of the population synthesis model 
such as /{, and (TO2) / (toi) infiuence the result at best by 
0(50%). These scaling relations allow us to estimate the 
merger rate implied by simulations {R) via a surrogate 
(i?) based on the number n of merger events seen in 
simulation: 



R = sn/T = 7.6 x 10 
~0.059yr-i 



M 



N {mi + TO2) 



(4) 



1 



N 1 + /b[(TO2) / (toi)] 



[Unless otherwise noted, "tilde" quantities refer to esti- 
mates for the corresponding "normal" quantity for one 
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individual simulation, based on only information from 
that one simulation.] 

Given unlimited computational resources, this simula- 
tion could be repeated many times to measure the merger 
rate R implied by this model to any accuracy desired, or 
equivalently measure the average number /i of of merger 
events expected from this simulation 

fj.= RT/s. (5) 

In practice, each simulation is performed once; there- 
fore, the simulated number of merger events (n) is only 
statistically correlated to the average number of events 
expected (/z), with relative probabilities p(n|/i) of any 
one simulation producing n events given by the Poisson 
distribution 

(6) 
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On average our estimates of the R — sn/T and equiva- 
lently of the number of events jl — n will agree with the 
true properties of the underlying simulation: averaging 

over many trials, (j^ = R and (/2) = /z. However, 

an estimate based on any one specific simulation should 
differ from the mean by a characteristic relative amount 



(logi? - logi?)2^ = V((logA - log^)') (7) 
~l/77Zln(10) ~ l/\AIln(10) 

2.2. Parameters varied in archives 

Our extensive experience in modeling of binary com- 
pact objects with StarTrack clearly indicates that there 
are seven parameters that st rongly influence c ompa ct 
object merger rates (see e.g., iBelczvnski et al.l ([2002:)): 
the supernova kick distribution (3 parameters cti.2 
and s describing a superposition of two independent 
maxwellians) , the massive stellar wind strength w (1), 
the common-envelope energy transfer efficiency a A (1), 
the fraction of mass accreted by the accretor in phases 
of non-conservative mass transfer fa (1), and the binary 
mass ratio distribution, as described by a negative power- 
law index r (1). We allow the dimensionless parameters 
a A, /a, w and s to run from to 1; the dimensionless r 
can be between and 3; and flnally we vary the disper- 
sion of either component of a bimodal Maxwellian ai , <72 
from to 1000 km/s. Additionally, motivated by recent 
observa tions suggesting pulsars with m asses near 2Mq 
(see,e.g. iNice et al.ll200a ISplaver et al.ll2002[ ) , we assume 
the maximum neutron star mass to be 2.5Mq'^. 

To improve our statistics, we combine results from sev- 
eral different databases of simulations. The most exten- 
sive database samples cti S [200, 1 200] and a-2 € fO, 200] 



very de nsely and was developed bvlO'Shaughnessv et alj 
(|2005aD and lO'Shaughnessv et all (|2005bl ) A second 



archive, significantly less dense due to computational re- 
source limitations, allows both dispersions to run uni- 
formly from to 1000 km/s. [Additional archives in- 
clude, for example, a set chosen to better-sample kick pa- 
rameters that best correspond to observations of pulsar 

^ Such a high neutron star mass converts many merging bi- 
naries we would otherwise interpret as BH-NS or even BH-BH 
binaries into merging NS-NS binaries, driving down the BH-NS 
and BH-BH rates significan tly from the distributions shown in 
lO'Shaughnessv et al.l l l2005bl ). as also described in §,4. 
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Fig. 1. — Scatter plot of the two bimodal kick velocity disper- 
sions (Ti , (72 for the population synthesis archives used to evaluate 
the merger rate of NS-NS binaries [denoted NSNS(m)]. The strong 
bias is introduced by incorporat ing high-density and computa tion- 
ally expensive simula t ions fr om lO'Shaughnessv et al.l (|2005al ) and 
lO'Shaughnessv et al.l ll2005bl '). 



rope r motions lArzoumanian et all (|1999f) : iHobbs et all 
2005( ).] Consequently, as shown in Figure[TJ our archived 

results do not uniformly sample the kick-related parame- 
ters through this range. This irregular sampling has two 
effects. First, having irregular sampling of the model pa- 
rameters effectively corresponds to imposing non-flat pri- 
ors on these parameters and hence biasing the resulting 
distribution function of merger rates that come directly 
from the database of runs; see Appendix \X\ and in partic- 
ular Figure IA6I Therefore it is even more important to 
develop the fits, which then allow us uniform sampling 
of the parameter space, and hence the derivation merger 
rate distributions assuming flat priors, as our intention 
is. Second, certain kick combinations are relatively un- 
dersampled, which likely plays a role in the relatively 
poor global convergence of fits for physical parameters, 
as described in § |31 Nonetheless, our sampling rather 
thoroughly explores the most p hysi cally likely regimes 
sugges ted bv lHobbs et al.l ()2005D and lArzoumanian et al.l 
(fT999l) . 



2.3. Event identification 

Most events of physical interest are uniquely and fairly 
unambiguously identifiable within the code. Type 11 and 
Ib/c supernovae events are distinguished by the presence 
or absence of a hydrogen-rich envelope at the supernova 
event. We also record DCOs which merge, so we can un- 
ambiguously determine the number of BH-BH, BH-NS, 
NS-NS, and WD-NS merger events which occur in a sim- 
ulation. [These event classes will be denoted BHBH(m), 
BHNS(m), NSNS(m), and WDNS(m) for brevity hence- 
forth.] Since the code also tracks binary eccentricity, for 
example, we can identify those WD-NS binaries which 
end their evolution with a non-zero eccentricy [denoted 
WDNS(e)]. 

When constructing our archived population synthe- 
sis results, we do not record detailed information about 
the nature and amount of any mass transfer onto the 
first-born NS. We therefore cannot determine the degree 
to which pulsars are recycled. However, we do record 
whether any mass transfer occurs. Thus for the purposes 
of identifying a class of potentially recycled ("visible") 
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wide NS-NS binaries [denoted NSNS(vw)], we assume 
any system undergoing non-CE mass transfer recycles 
its NS primary. 

2.4. Practical Archive Generation and Resolution, with 
Partitions and Heterogeneous Targets 

The accuracy to which each population synthesis event 
rate prediction is known, ~ 1/ y/n, is uniquely set by the 
number of events n seen within that simulation. For this 
reason, lO'Shau ghncssv et al.l ((2005a) (i) designed their 
population synthesis runs to continue until a fixed num- 
ber of events were seen and (ii) used results only from 
such targeted simulations, where a minimum number of 
events was guaranteed. In contrast in this study, given 
limited computational resources and a wide range of tar- 
gets for which predictions are needed, we extract all pos- 
sible information from each simulation: whenever possi- 
ble, we make an estimate of each event rate of interest. 

However, it is important to note that most of our sim- 
ulations employ some degree of accelerating simplifica- 
tion which can bias estimates. To give the most ex- 
treme example: the most accurate estimates for the BH- 
BH merger rate come from population synthesis runs 
which evolve only a subset of possible progenitor bina- 
ries by using partitions. This subset has been shown 
to include the progenitors of the vast majority of BH- 
BH binaries, but very few progenitors of NS-NS binaries 
and other less massive DCO s (for more information, see 
lO'Shaughnessv et al.ll2005ai ; we continue to employ the 
same partition they devised). Similarly, the vast major- 
ity of simulations used to study NS-NS event rates (i) use 
a similar partition to reduce contamination from white 
dwarf binaries and (ii) terminate the evolution of any 
binary immediately after any WD forms. These strong 
biases make data from these two types of simulations in- 
appropriate for use in, for example, estimating the WD- 
NS merger rate. For this reason, the number of popula- 
tion synthesis archives A'' available to make predictions 
varies significantly across the various target event types; 
see Tabled 

Furthermore, the accuracy each simulation can provide 
varies considerably: unlike earlier studies, the number of 
event samples n in each archive is not guaranteed to be 
greater than a minimum value. However, usually we have 
at least one event for each event type: Ns+, the num- 
ber of unbiased population synthesis results containing 
one or more events is usually very close to TV^, the to- 
tal number of unbiased population synthesis simulations 
available. And in most cases a significant proportion of 
our sample contains enough events to determine the rate 
to better than 30% (i.e., n > 10). 

3. FITTING: PRINCIPLES AND TESTS 

Given the prohibitive computational demands of di- 
rect population synthesis simulations described in § [21 we 
use fits based on archived results of population synthe- 
sis runs as a surrogate for repeated detailed simulations. 
Confidence in our results is therefore tied intimately to 
confidence in the quality of these fits. However, even low 
order fits in seven dimensions involve many parameters: 
to fit any nontrivial function, we must fit roughly a hand- 
ful of data points per parameter. To build confidence, we 
must show that the fit order chosen adequately describes 
the data without overfitting. More delicately, since these 



fits are used in the derivation of empirical constraints on 
the population synthesis mod els in our companion paper 
()0'Shaughnessv et al.|[2007br ). we must also be able to 
show that the key end product, the "constraint-satisfying 
model region" does not depend sensitively on the fit de- 
tails or on random accidents in the data (i.e., any dif- 
ferent Monte Carlo realization of our simulations should 
yield the same result). Finally, we note th at our com- 
panion paper (lO'Shaughnessv et al.l[2n07bl ) wiU impose 
four independent constraints. In order to have more than 
90% confidence that all models that satisfy all constraints 
are indeed inside the intersection of the four constraint- 
satisfying regions, we must at a minimum show that each 
individual constraint-satisfying prediction contains more 
than 0.90^/^ ~ 0.974 of the models that truly satisfy that 
constraint. 

3.1. One-Dimensional Model 

Though we intend to build confidence explicitly in 
our ability to make predictions on the basis of seven- 
dimensional fits, some of the required principles, no- 
tation, and tools are best illustrated through a one- 
dimensional example. Thus in this section we explore 
an arbitrarily-chosen but known "population synthesis" 
model. Depending on one real parameter x with < x < 
1, this model on average will produce ijl{x) merger events 
out oi N = 10^ binaries, with /x(x) chosen arbitrarily as 

^(:e) = 10^^'") = 5 + 20x2 . (8) 

Specifically, in this model we assume: all stars are 
born in binaries (/f, = 1); that the smaller companion 
star has mass randomly distributed from the hydrogen 
burning limit to the primary's mass, corresponding to 
{^2) I {1^1) = 0.5; and that star formation extends over 
an interval T = lOGyr. Under these specific circum- 
stances, as discussed in fj2] [Eqs. (|5|3p ]. a model which 
produces fJ.{x) events on average out of 10'^ binaries cor- 
responds to a merger rate 

R{x)^s^i{x)/T (9) 

~4 X 10"V(a;)yi-"^ • (10) 

An individual simulation of this model will produce some 
number n of events, with n statistically related to /i(a;) 
via the Poisson distribution [Eq. ([6])]. 

To estimate the merger rate as a function of x, we 
perform several simulations [Ns = 11) using different 
parameter values 1)/10 for index val- 

ues a = 1 ... 11, obtaining results Ua and thus estimates 
Ra = sria/T as shown in Figure [21 To better fit the 
merger rate over the many orders of magnitude that ap- 
pear in practice, we choose our fit R{x) = io*^("^)s/r 
to minimize the logarithmic difference between it and 
simulations, weighted by the relative statistical uncer- 
tainties of each simulation, l/n(lnl0)2 from Eq. (O. 
[Here and henceforth we use "hat" symbols to denote 
our fit quantities; unlike "tilde" quantities, which esti- 
mate properties of an individual simulation, the "hatted" 
quantities generally depend on the results of all simu- 
lations simultaneously.] Specifically, to generate a g-th 
order polynomial fit log Rq {x) for the logarithm of the 
merger rate, we choose the coefficients of the polynomial 
Mq{x) = log Rq{x)T / s to minimize the "average square 
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deviation divided by characterist deviation, per degree 
of freedom" 
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where N<q is the number of (/th-order basis polynomials 
(g + 1, in one dimension). The fit which minimizes Xg is 
the (gaussian) maximum- likelihood estimator for R{x), 
the best single estimate possible.^ Figure shows the 
results of this minimization for the linear (g = 1) and 
quadratic (q — 2) fits: 

log i?i=log(s/r) + 0.49 + 0.93X (12) 

logi?2 =log(s/r) + 0.35 + 1.49a; - 0.547x^ (13) 

Both performing simulations and fitting introduce er- 
rors, which we can quantify by the mean square deviation 
between the exact merger rate logi? and its fit log-R^: 

-,1/2 



In 



dx\\og{Ry{x)/R{x))\' 



(14) 



For example, the value of /i is 0.075. Roughly speaking, 
Iq estimates the rms uncertainty in the fit: Rq should lie 
within a factor ~ 10^^' of R{x). For example, Ri should 
be accurate to within a factor roughly 10='=0-'''^^ ~ 1.2=*=^, 
as confirmed by Fig. Similarly, the values of x that 
correspond to a given merger rate should be uncertain to 
roughly O {I / [dlogR/ dx)) ~ 0.1. 

As noted previously, an extremely high level of con- 
fidence is needed in any prediction of a constraint- 
satisfying region. Let us use V{G, /) to denote the set 
of points that / maps to an interval G = {(/mim 5max}i 
so V{G,f) — {x\f{x) £ G}; and C is an observation 
we demand our simulation reproduce; then the fraction 
of V{G, log R) - the population synthesis models which 
truly satisfy the constraint - which are inside V{G, log R) 
- the set of population synthesis models which are naively 
predicted to satisfy the constraint, based on the fit and 
the original constraint - is given by r+(C, log R\C, log R) 
where is generally defined by 



r+{A, log R\B, log R) 



\V{AlogR)nV{BlogR)\ 
l^(^logi?)| 



(15) 

for arbitrary intervals A, B; here \ V\ is the volume of V. 
In the example shown in Figure [2l the fraction r_|_ of the 
truly constraint-satisfying area that is predicted to sat- 
isfy the constraints is only 78%; if this constraint had to 
be combined with four other constraints of similar qual- 
ity, fewer than (78%)'* ~ 37% of all constraint satisfying 
points would lie inside the intersection of all four naive 
predictions. 

To increase the fraction of constraint-satisfying points 
included in a prediction - in the notation of this pa- 
per, to find a volume V which contains the "naive" 

^ We have also appled a maximum-likelihood estimate based 
on Poisson rather than (approximate) gaussian errors. For the 
problems explored here, for which many merger events are usually 
available, we see no significant difference between the results of this 
more statistically appropriate method and the conventional and 
pedagogically far simpler gaussian maximum-likelihood method. 

Of course, regions where the fits perform significantly worse 
can exist; here, the fit differs from logi?(a;) by ~ 6.5/i near a: = 0. 
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Fig. 2. — Illustration of how small errors in fits can produce 
large uncertainties in the constraint-satisfying region, using a one- 
dimensional model. Solid line: the "true" function log R{x) = 
log fi{x)s/T [Eq. its}], expressing the log of the average number 
of events expected for any x. Points: results of "experiments" 
for eleven values of the parameter x. Dotted lines: linear and 
quadratic least-squares fits to the data. Dark shaded boxes: (y- 
axis) the "constraint" interval C = [—5.8, —5.4]; (x-axis) the range 
of X which truly satisfies the constraints [i.e., all x so logi?(x) G C; 
this region is denoted V {C,log R) in the notation of this paper]. 
Light shaded box: the interval predicted to satisfy the constraints, 
based on a linear fit [denoted V{C, logRi) in the notation of this 
paper] . 



prediction V{G, log Rq) C V but contains more of the 
truly constraint-satisfying points ^(C, logi?) - the sim- 
plest option is simply to increase the size of the "con- 
straint" interval. A larger constraint interval C* con- 
taining the original would automatically include more 
points in its predicted region (now V{C* ,logR)). While 
this prediction loses some of its reliability (i.e., we can- 
not ensure that only constraint-satisfying points lie in 
this region), we increase its robustness: by a good 
choice of G* we can almost guarantee that all constraint- 
satisfying points are included. Specifically, if we change 
the interval C = {cmin, Cmax} to the broader interval 
C*{I) = {cmin — i, Cmax + /} ^ that is, if wc incrcasc 
the constraint region by the size of the characteristic er- 
ror in the fit - then any point x which truly satisfies 
the constraints (in our notation, x G V{C\log R)) al- 
most certainly lies within the larger "predicted" volume 
V{C* ,logR). Formally, assuming the error at any point 
|(51ogi?| is likely less than /, we conclude that logR{x) = 
Slog R{x)+ log R{x) > S log R{x) + Cmin > ~I + Cmin and 
similarly that logR{x) ^ I + Cmax- In the case shown in 
FigureHl this error-widened prediction V{G*{Ii),logRi) 
includes both a higher fraction of constraint satisfying 
points (r_|_ ~ 92%) and a higher fraction of points which 
are not constraint satisfying (20%). 

All these calculations, however, depend on the fit or- 
der: higher fit orders q often allow more degrees of free- 
dom with which to reproduce the true function log R{x) 
and therefore estimate y(C, logi?). For example, as- 
suming arbitrarily small (constant) errors and arbitrar- 
ily dense sampling Xa, the best possible g-th order 
polynomial least-squares fit is the orthogonal projection 
Pq log R of log i? onto the space of g-th order polynomi- 




Fig. 3. — Illustration of the balance between truncation and 
sampling error in a one-dimensional fit. 



als, given in one dimension by 

9 pi 

{Pg\ogR){x)=Y,Mx) / dyMy)^ogR{y) 
1=0 -^0 

for 4>i{x) any set of orthonormal basis polynomials of 
degree < I (e.g., (j)i{x) = \J21 + lP/(2x — 1) where 
Pi{x) are Legendre polynomials). Therefore, the mini- 
mum possible error of a gth-order polynomial fit is the 
magnitude ||Pj_qlogi?|| of the difference Pj_qlogi? = 
logi? — Pq log R between logi? and this best fit, where 

II/IP = Jq dx\f{x)\'^. As the fit order increases, this 
"truncation error" 



et,9 = \\P±q^ogR\ 



(16) 



decreases dramatically, as shown in Figure [3l However, 
the total error Ig in the fit reflects a balance between (i) 
decreasing the "truncation error" with more degrees of 
freeom and (ii) increasing the "sampling error" as fewer 
data points are available per degree of freedom. For this 
reason, an optimal fit order exists, which provides the 
smallest error in R{x), indicated by a minimum in Iq. 
This best fit order will provide the most reliable estimate 
oiV{C, logi?). 

In principle, this one dimensional example contains 
all the core concepts we will apply to fitting seven- 
dimensional population synthesis results: the sampling 
uncertainty in each sample Ua of "simulated merger 
events" ; the fit and its inaccuracies; the need to increase 
the size of the constraint interval and thus the "pre- 
dicted" constraint-satisfying region to include a larger 
fraction of constraint-satisfying events; and a specific 
procedure for doing so, based on the characteristic errors 
of the fit. In practice, however, we must substantially 
flesh out this outline in two ways: (i) we must add a re- 
liable way for estimating the fit error and selecting the 
optimal fit order without using /, since its definition in 
Eq. (fH)) uses the fit itself and requires exact, a priori 
knowledge of R{x) (the function being fitted) that is not 
available in practice; (ii) we must extend the concepts de- 
veloped here to seven dimensions, given the dependence 
of our population simulations. 

3.2. Seven-Dimensional Analog 



By virtue of existing in seven dimensions, the full popu- 
lation synthesis fitting problem is qualitatively different 
than the one-dimensional problem even with the same 
number of points and characteristic error. In the one- 
dimensional case, the many densely packed and relatively 
accurate samples of R permit very accurate fits to log R; 
barring rapid variation in logi? with x, one-dimensional 
fits can easily have accuracies considerably exceeding the 
limiting accuracy of any individual measurement. On 
the other hand, in seven dimensions the same number 
of data points are much more loosely spaced (geometri- 
cally, the characteristic spacing is proportional to N~'^/'^) 
and must constrain many more degrees of freedom at 
the same polynomial order, yielding much less accurate 
fits. For this reason, rather than attempt to rescale the 
number and uncertainties of our population synthesis 
simulations to generate a tractable one-dimensional ana- 
log, we instead demonstrate convergence using a seven- 
dimensional model that best captures the relevant fea- 
tures of our population synthesis data. 

Explicitly, in this seven-dimensional toy model (i) the 
logarithm of the simulation size logiV is gaussian dis- 
tributed around —4.5 with an order of magnitude stan- 
dard deviation, omitting simulations smaller than 10"^ 
and larger than 10^; (ii) each set of simulation param- 
eters Xa is chosen by uniformly selecting seven random 



numbers < 



^fc=1...7 



< 1; (iii) the number of mergers 



Ua observed in any particular simulation parameterized 
by the seven parameters Xa is drawn from a Poisson dis- 
tribution with mean ^{xa) — p{xa)N, where the mean 
number of mergers per binary p{x) is defined by 

pix) = fl{x) /N = ^ (17) 

and thus (iv) implies a merger rate 

i?(x) = p(x)(0.04yr-i) (18) 

(again assuming /{, = 1 and (7712) / (mi) = 0.5). More 
specifically still, to compare with the array of 0(2000) 
simulations of merging NS-NS binaries, we generate a 
"model archive" of Ns = 2000 randomly realized analogs, 
each defined by their parameter choices Xa , specific simu- 
lation sizes Na , observed merger numbers , and merger 

rates Ra = OMyr-^ina/Na). 

This seven-dimensional model qualitatively reproduces 
the distribution of merging double neutron star [NS- 
NS(m)] simulations in {n, N) and i? [Fig. [4]. Quanti- 
tatively, this toy model is on average as accurate as our 
population synthesis simulations: the root-mean-square 
relative sampling error as (the "expected" standard de- 
viation in the data, given the limiting error produced by 
sampling at each point) 



1 ^= 



(19) 



of our population synthesis simulations o^e = 0.25 [see 
Table [2] agrees with the corresponding value as = 0.24 
[see Table [T] for this seven-dimensional analog. Fur- 
thermore, in both cases this average accuracy ctb is still 
slightly smaller than the range in log i?, as characterized 
by (Jdd (the "data-data" standard deviation, measuring 
the range of the distribution): 

"2 - ^ E[(i°g^")' 



'DD 



logi?\' 



(20) 
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Fig. 4. — Comparison between our toy model (left) and sim- 
ulations of merging NS-NS binaries (right). Top panels: scatter 
plot of the simulation size A'^ versus the number of merger events 
n seen. Bottom panels: distribution of merger rates seen in simu- 
lation. In the right panels, no corrections are being made for the 
highly irregular patterns in which N and n were chosen, and which 
introduces unavoidable biases in the distribution of log(R,). 



{o'DD = 0.73 for real simulations; ano = 0.42 for our 
toy model) Finally, as can be seen by comparing A^^, 
the number of population synthesis simulations, to N^, 
the number of simulations with at least one event [Table 
[1], a significant fraction 0(10% — 50%) of population 
synthesis simulations are totally unresolved; our seven- 
dimensional toy model has a similar fraction (20%) of 
unresolved simulations. 

3.3. Error Estimates and Convergence 

Using this concrete and highly realistic seven- 
dimensional model problem described in ^S.2\ we can not 
only conclusively demonstrate that fairly accurate fits are 
possible for population syntheses, but we can also find a 
reliable estimate of the rms error in the fit, an estimate 
that can be reliably applied to surround all constraint- 
satisfying points inside a "predicted" volume. Precisely 
as in the one-dimensional case, given a fixed archive 
of seven-dimensional simulations, we can estimate log R 
with logi?^, the unique polynomial of order < q that 
minimizes Xq^ using any integer q [Table [Tj. With more 
degrees of freedom, higher-order fits can much more eas- 
ily reproduce the data, so their mean square difference 
from the data ctdp 

~ 1 2 



2 



\ogR{Xa) - logi?a 



(21) 



(a "data-fit" comparison) decreases monotonically; see 
the fourth column of Table [1] But as these high order 
polynomials increasingly fit every sampling-induced flu- 
cuation in the data, high-order fits increasingly deviate 
from the exact solution: the exact rms error Iq is also 
shown in Table [1] [calculated by comparing the fit to the 
known seven-dimensional model R{x) from Eqs. (??)], 
showing Iq has a minimum near an optimal fit order. 
Near this optimal fit order, the characteristic fit error Iq 
is substantially less than the characteristic range of the 
data {auo) "in fact, even less than the average statisti- 
cal error of the sample ge- In other words, though this 
fit is accurate to only 20% — 30% at any point (based 
on 10^' - 1 for q = 1. . .4), this accuracy is sufficient to 



Table 1. Behavior of Sample Fit as Order Increases 
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J 


/ 


O'DF 


O'DD 


0"B 





60.7 


0. 


0.462 


0.421 


0.421 


0.241 


1 


4.16 


0.373 


0.138 


0.228 


0.421 


0.241 


2 


1.3 


0.104 


0.0813 


0.203 


0.421 


0.241 


3 


1.61 


0.0475 


0.0929 


0.198 


0.421 


0.241 


4 


2.4 


0.105 


0.143 


0.182 


0.421 


0.241 



apply constraints, given the many orders of magnitude 
range spanned by R{x). 

But while accurate fits exist, the tools presented hence- 
forth provide few ways to identify them and particularly 
estimate their error without resorting to knowledge of 
the function being fit. For example, as we increase the 
available number of degrees of freedom, the root-mean- 
square difference aoF between the fit and the data de- 
creases monotonically to zero. Even at or near the op- 
timal fit order, in our experience both ctdf and ap, the 
expected sampling-induced root-mean-square difference 
between data and predictions [Eq. (fT9|) ]. correlate only 
weakly with the true error in the fit {Iq). A good fit can 
be identified^ using the weighted squared difference Xq^ 
which does have a minimum value for a good fit.^ How- 
ever, to understand and estimate the fit-induced error we 
must compare fits of successive orders with 



Jq^ 



lOg{Rq/Rq-l) 



1/2 



(22) 



This comparison between fits of different orders arises 
naturally from understanding how the error in the fit 
Iq reflects a balance between (i) an increasing ability 
to fit logi? exactly with more basis polynomials and 
ever smaller "truncation error" , decreasing the magni- 
tude I \P±q log R\ I for P±q the and (ii) a decreasing statis- 
tical accuracy with fewer independent degrees of freedom 
Ns — N<q available to constrain the N<q independent co- 
efficients in the expansion. These two terms contribute 
independently to the error at any fit order, allowing us to 
estimate Iq by an incoherent superposition of "truncation 
error" and statistical error without explicitly calculating 
the best fit log Rq itself: 



\P±q 10gi?| 



4 



N. 



<q 



- N. 



<q 



with 



N. 



<q- 



{d + q)\ 
d\q\ 



(23) 



(24) 



Here c? = 7 is the dimension in which our functions are 
defined and Pig projects functions perpendicular to the 
space of basis polynomials of order less than or equal to 

^ The best way to identify a good fit, given sampled data, is 
with a blind test, where data points not involved in the fitting 
process are compared against the fit. We have explored using blind 
tests points and calculating 51 the analagous weighted squared 
diff'erence between the fit and these test points. In our test cases, 
J< 

blind tost. 
^ In fact 

the sampling statistics allow, when x 



Xq also allows us to identify when the fit is as good as 



1 
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q [cf. Eq. p6)) ]: for example, if logi? is a polynomial, 
P±q selects those coefficients of log R of order less than q. 
In this and other cases, this expression very accurately 
reproduces Ig ; see for example Table [1] 

When sampling errors are relatively small (i.e., Ua ^ 
1, so cte ^ 0), each successive approximation order 
log Rq ~ Pq log R is found by orthogonal projection of 
log R; therefore, for low fit orders the truncation error is 
just the magnitude of the next-highest-order correction 
to be applied, Jq-. 

||Pi,l0gi?|| ~ ||l0gi?,+ l -l0gi?,|| = Jq+l 

In other words, for low q we expect Iq ~ Jq+i, as con- 
firmed in Table [TJ 

From our experience with this and other model prob- 
lems, the true rms fit error (Iq) is quite generally close to 
the change in the fit to the next order (J^+i) and from 
the previous order (Jq). In particular, the difference be- 
tween the current and next-lowest-order fit (Jq) has a 
minimum similar in magnitude and location to the mini- 
mum of Iq; see for example Table [TJ We therefore use the 
order at which Jq is smallest to select the optimal fit or- 
der and its value as the characteristic root-mean-square 
error in the fit. 

3.4. Volume estimation 

This error estimate provides the key tool needed to reli- 
ably and algorithmically surround almost all of the points 
with merger rates R consistent with constraints, using in- 
formation only about the fit, to the high level of accuracy 
needed to trust the results of multiple constraints. For 
example, observations of Galactic binary pulsars suggest 
double neutron star binaries merge in the Milky Way 
merge at a rate between Cmin = log(3 x 10~^yr~^) and 
Cmax — log(23 X lO^^yr^^). This constraint range turns 
out to lie on the high end t ail of what our model simu- 
lations produce ()Kim et al.l 2003; 0'Shaughncss v~et al.l 
I2007bfl . If our seven-dimensional model exactly de- 
scribed double neutron star merger rates, then only a 
small fraction of model parameters x £ V{C,logR) (5% 
by volume of the unit seven-dimensional cube, corre- 
sponding to \x\ < 1.04) reproduce observations. How- 
ever, even though these merger rates correspond to the 
largest possible and therefore best-sampled predictions 
from population synthesis, the fit remains sufficiently in- 
accurate so that only 1 — r_|_ ~ 18% of all constraint- 
satisfying points would not be included in a prediction 
y(C, logi?) based on the fit alone. But based on Table 
[U we can estimate the characteristic error in our fit by 
J ~ 0.05. Therefore, when we compensate for this un- 
certainty with a wider constraint interval C* = C*(J), 
we much more accurately bound the set of constraint- 
satisfying points: 1 — r+(C*, log -R) = 5%. 

Strictly, the fraction of constraint-satisfying points 
that will be encompassed with a larger constraint inter- 
val C depends strongly on the width and placement of 
the original interval C - intuitively, an extremely narrow 
constraint requires extremely high-precision reconstruc- 
tion of i?(a;), in turn possible only at high merger rates R, 
where many events should be seen in simulation. How- 
ever, the case shown here, with a very tight constraint 
(with only 5% by volume of parameters satisfying the 
constraint), is substantially more difficult to satisfy than 



most individual observational constraints that can be ap- 
plied in practice. Comparing our experience with model 
problems to the weak constraints available, we are confi- 
dent that fewer than 5% of constraint-satisfying systems 
will be omitted from V{C, log R) when these volumes are 
constructed on the basis of real population synthesis sim- 
ulations and corresponding observations. [The results of 
our broader array of model problems are available online 
in .O'Sh aughncssv ct al. (2007a).] 

4. FITTING: POPULATION SYNTHESIS 
PREDICTIONS 

Almost exactly the same fitting techniques presented 
above are applied to our population synthesis data: we 
perform a weighted least-squares fit of polynomial-like 
basis functions over our seven-dimensional space, for 
each of the event rates of interest [BHBH(m), BHNS(m), 
NSNS(m), NSNS(vw), WDNS(e), WDNS(m), SNIb/c, 
and SNII, as introduced in ^2.3| . For each fit, we eval- 
uate the fit quality of several different polynomial or- 
ders to our data. To minimize the possibility of using 
more parameters than allowed by the number of our data 
points, we choose as "best" fit that order which mini- 
mizes the relative difference between fit orders Jq (see 
below). However, we know that the merger rate R{x) 
must satisfy certain symmetries, based on the manner in 
which we represent the kick velocity dispersion as a sum 
of two maxwellians (i.e., we can switch labels associated 
with the two distributions); therefore we limit the man- 
ner in which those parameters associated with kicks can 
enter into the distribution (see the Appendix). 

Table [5] summarizes the properties of the least-squares 
fits we applied to our archived population synthesis re- 
sults. The first column provides a brief label for the 
fit, as described in greater detail in the text. The next 
two columns summarize the amount of available informa- 
tion contained in our population synthesis archive: N is 
the number of population synthesis models with unbiased 
data (i.e., where all plausible progenitors for the target 
event have been included), whereas is a smaller num- 
ber of models with unbiased data containing one or more 
events (i.e., for which an estimate of the rate, rather than 
merely an upper bound, is possible). The next block of 
five columns describes properties and diagnostics of a 
weighted least-squares fit applied to our data. The first 
two columns, q and N<q, merely indicate the polynomial 
order and number of degrees of freedom involved in our 
fit {N<q differs from N<q introduced earlier due to the 
symmetry requirement described above) . In all the cases 
shown here, the optimal polynomial order produces far 
fewer degrees of freedom than population synthesis sim- 
ulations (i.e., N<q <C Ns^+). 

The next three columns provide critical diagnostics of 
our fit: x^ CTDF, and Jq [Eqs. (|11I19I22|) ]. All three 
columns roughly measure our "goodness-of-fit" ; based 
on our experience with model problems and our un- 
derstanding of the values these quantities should take 
for fits dominated by statistical and truncation errors, 
all three universally indicate that truncation error dom- 
inates - that is, that our low-order polynomial basis 
functions are not sufficiently general models to match 
the exact merger rate R{x) implied by our population 
synthesis simulations. For example, the limiting value 
of Jq is significantly above the level of sampling error 
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Table 2. Statistics of population synthesis simulations for event rates 



Type 


Ns 






q 






J, 


O-E 


O'DD 


77(10) 


BH-BH(m) 


2930 


1201 


141 


3 


15.5 


0.77 


0.67 


0.25 


0.78 


12% 


BH-NS(m) 


2533 


1334 


141 


3 


11.1 


0.53 


0.40 


0.23 


0.71 


2% 


NS-NS(m) 


2803 


2382 


141 


3 


18.4 


0.37 


0.22 


0.14 


0.63 


5% 


NS-NS{vw) 


1325 


1087 


141 


3 


16.7 


0.39 


0.37 


0.13 


0.76 


8% 


WD-NS{c) 


1770 


1564 


141 


3 


12.5 


0.34 


0.20 


0.14 


0.57 


10% 


WD-NS{m) 


1770 


1658 


141 


3 


16.4 


0.34 


0.19 


0.13 


0.45 


11% 


SN Ib/c 


1482 


1482 


141 


3 


7.8 


0.07 


0.06 


0.02 


0.11 


n/a 


SN II 


1482 


1482 


141 


3 


6.0 


0.04 


0.02 


0.03 


0.17 


n/a 



<^e/ \/Ns^+/N<,q suggested by the second term in Eq. 

- that is, above the level of error expected when av- 
eraging Ns^+/N<q simulation samples per degree of free- 
dom, each with characteristic error ~ as- This high level 
of error suggests fit accuracy is limited by truncation er- 
ror - inability to fit the exact form of R{x) with our 
basis polynomials - rather than sampling uncertainties 
at each point. ^ Additionally, as illustrated in our dis- 
cussion of seven-dimensional model problems, the latter 
two columns provide the logarithmic uncertainty in our 
fits; for example, fits with J ~ 0.33 are known to within 
a factor two with one-sigma confidence. These columns 
should be contrasted with the next two, which provide 
the average sampling-induced uncertainty in the data 
[Eq. p9)) ] and the characteristic range of merger rates 
ajjD seen in simulation [Eq. ([20])]. Our fits remain use- 
ful so long as their uncertainties are much less than the 
range of logi? (i.e., so long as J < aDo)- As seen in Ta- 
ble [21 our best fit satisfies this requirement (J < (Jdd) 
for each populations of interest. 

The last column provides the fraction rj of simulations 
with no events when more than 10 would be expected 
based on the fit. Based on an average uncertainty of 
a characteristic factor of two in the fit (our best fits 
have J ~ 0.33, corresponding to a merger rate known 
to within a factor 2 ~ 10"'), we would estimate that the 
fraction of simulations rj which produce no merger events 
when more than 10 should be seen based on our uncer- 
tain estimate of the event rate should be comparable to 
or smaller than the fraction of simulations which should 
produce 10/2 = 5 events but in fact produce zero, namely 
T] < p(0, 10/2) ~ 0.7% [Eq. dH)]. However, regions with 
the lowest merger rates will be undersampled and there- 
fore have characteristic uncertainties marginally larger 
(e.g., up to a factor 3 to 4), leading to ry of a few to ten 
percent. 

4.1. Results 

Supenovae: Given our choice of stellar mass interval 
probed in our simulation (i.e., mi > 4Mq), supernovae 
occur extremely frequently, providing us with superb 
statistics at low cost. However, our limited set of basis 
functions can only with difficulty reproduce the observed 
variation in supernovae rates: even though SN rates for 
models in our archive are at times known to 1% (i.e.. 

Additionally, when fitting supernovae rates as a function of 
population synthesis parameters for single stars, where only one 
of our parameters (wind strength) enters, we find the rate has 
a moderately complex functional form that requires high-degree 
polynomials to fit. We expect similarly complex behavior in the 
multidimensional case, and interpret the large correspondingly. 



involving 10"* or more sampled events), or 0.004 in the 
log, our optimal fit differs significantly from the data, by 
CT ~ 0(0.04) in the log. 

Nonetheless, as d iscusse d in the forthcoming 
lO'Shaughnessv et all (|2007bD . the supernovae rate 
remains a striking success of the StarTrack population 
synthesis code and our normalization conventions (e.g., 
M ~ 3.5Moyr~^). No matter what combination of pop- 
ulation synthesis parameters we choose, the predicted 
SN rates lie well within the observational constraints 
found bv lCappellaro et al.l |l999). 

WD-NS binaries: As with supernovae, white dwarf- 
neutron star binaries occur fairly frequently,allowing us 
to accumulate fairly good statistics over a broad range 
of population synthesis parameters. Additionally, based 
on the distribution of Ns versus n (i.e., as in Figure |A5|) . 
our sample shows no signs of systematic incompleteness: 
we appear to have covered the full range. Though our 
polynomial fits continue to introduce systematic error, 
the resulting fit behaves well throughout the range. 
BH-BH binaries: Double black hole binaries, in con- 
trast, occur extremely infrequently, especially with an 
assumed maxi mum NS mass of instead of 2. OA/©, 

as assumed in lO'Shaughnessv et al.l (2005a, b). [The in- 
efficient formation of coalescing BH-B H binaries is dis- 
cussed and explained in more detail bv iBelczvnski et "ahl 
(l2006bD.] N onethe less, by using s pecial-purpose parti- 
tions. lO Sha ughne ssv et al.l (|2005aD accumulated a large 
(~ 500 simulations) sample with good {n ~ 10) statis- 
tics on a large subset of parameter space; though these 
simulations assumed a maximum NS mass of 2Mq, we 
post facto changed the maximum NS mass to 2.5Mq. 
By adjoining the results of general-purpose simulations 
not assured of good statistics, this sample has since been 
enlarged by a factor ~ 5, with emp hasis on the same sub- 
set of p arameter space presented in lO'Shaughnessv et al.l 
(|2005a| ) (see Figure [ij . Fin ally, though the distribution 
of n versus Ns (Figure lAsp suggests the lowest BH-BH 
merger rates may not be very well-resolved, we have no 
reason to suspect we have any significant under-resolved 
region: simulations larger than 10^ binaries consistently 
produce several merger events (n ^ 0). 

Nonetheless, the BHBH(ni) fit is comparatively poor: 
even though the average simulation with binary black 
holes produces enough to fairly accurately determine the 
rate (i.e., ue — 0.25), our best fit is barely more accurate 
than approximating the average BH-BH merger rate seen 
in simulations (i.e., compare the characteristic fit error 
Jq ~ 0.67 with the range of BH-BH merger rates seen in 
simulation, aoD — 0.78). 
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Results: NS-NS and BH-NS binaries: Despite (or per- 
haps because of) a concerted effort to accumulate good 
statistics targeted specifically to these classes, fits for 
NS-NS and BH-NS event rates are statisticalUy implau- 
sible as measured by x^- Judging from Figure IA5[ 
the NSNS(m) and BH-NS (m) are comparatively well- 
resolved: longer simulations fairly consistently have a 
lower chance of producing n = results. For this reason, 
we strongly suspect some feature of these underlying rate 
functions are poorly-described by our basis functions; we 
intend to more thoroughly test this hypothesis (with bet- 
ter statistics) by comparing these fits to nonparametric 
estimates in a future paper. 

5. CONCLUSIONS 

To develop a more comprehensive understanding of 
population synthesis predictions and to allow those the- 
oretical predictions to be systematically compared with 
observations of the end products of high-mass single and 
binary stellar evolution, we have fit eight predictions 
from the StarTrack code over seven of its input param- 
eters. These fits are available on requests sent to the 
first author. In a companion paper. lO'Shaughnessv et al.l 
(|2007bl apply these fits along with estimates of their sys- 



tematic errors to discover robust constraints on the seven 
parameters that enter into population synthesis. Addi- 
tionally, we have demonstrated that in analagous model 
problems the constraint-satisfying region defined by us- 
ing these fits can, under appropriate conditions, very 
accurately trace the underlying constraint-satisfying re- 
gion. Finally, we have presented a thorough diagnostic 
formalism, including a large list of diagnostic quantities 
and tests (I , J,(Te,df.,dd,X^ iVi'^+i etc.), which can be 
applied to studies of fits to large archives of any (indi- 
vidual) population synthesis code. 
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APPENDIX 



SAMPLING, PRIOR DISTRIBUTIONS AND THE LIKELIHOOD OF RESULTS 

In Figure lA6l we compare the distribution of merger rates derived directly from our current database of simulations 
(dotted lines) to the distribution of merger rates derived from the corr e spondi ng fit (solid lines) obtained using the 
methodology described here. Unlike Figure 3 of lO'Shaughnessv et al.l (|2005bl ). in which the simulations being fit 
had parameters drawn randomly from the entire range of parameters allowed, in this study the simulations were not 
randomly distributed through the entire parameter range; see for example Figure [TJ In effect, the rate distributions 
drawn from simulations and fits assume significantly different priors for population synthesis model parameters, the 
former favoring lower kicks. Therefore, although the distributions appear different in shape, they refiect the same 
physical process, just with different priors regarding what model parameters are likely. In particular, the fitted rate 
results represent the unbiased distributions of rates with uniform coverage of the seven-dimensional model parameter 
space. Furthermore these distributions fully account for uncertainties in the fits, as described in 21 
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Fig. A5. — Log-log scatter plot of the number of events seen in a sample against the population synthesis sample size. Also shown is 
a diagonal line corresponding to the 99% lower bound on the fit-predicted rate distribution, translated from physical rate into expected 
number of events seen per unit population synthesis sample size. For completeness, models with zero events are shown with logj^Q n as — 1. 



log(R) log(R) log(R) 

-8. -7. -6. -5. -4. -3. -8. -7. -6. -5. -4. -3. -8. -7. -6. -5. -4. -3. 




-7. -6. -5. -4. -3. -2. -7. -6. -5. -4. -3. -2. -7. -6. -5. -4. -3. -2. 
log(R) log(R) log(R) 



Fig. A6. — Probability distributions for formation rates for different types of binaries in the Milky Way. Dotted lines: the distributions 
of rates derived directly from our database of simulations, which does not regularly sample the model parameter space; cf. Figure [T] Thin 
solid curves: the distribution of different formation rates derived using fits, using precisely the same sampling as the raw data (dotted line). 
Heavy solid curve: the distributions derived from our fits to the simulations when the model parameter space is sampled uniformly. In all 
cases, the formation rate data is smoothed with a gaussian of standard deviation 0.1. It is evident that: (i) the quality of the fits is rather 
good; and (ii) the change in the shape of the distributions is primarily because of the different priors between the dotted and thick-solid 
lines. Last we note that: (i) the rate results for binaries involving BHs are most uncertain because of the assumed maximum NS mass of 
2.5Mq that reduces the overall formation rates for BH binaries (compared to an assumed maximum mass of 2Mq); (ii) all rates quoted in 
this studies are for Milky- Way like galaxies and we do not account for any differences associated with elliptical galaxies. 



To phrase the same statement more abstractly, our prior expectations about the relative likelihood of population 
synthesis parameters (e.g., kicks) influences our expectations for the relative likelihood of different merger rates. If 
at any point our knowledge of binary properties, evolution, and compact object kicks improves to a degree such that 
we can confidently move away from flat priors into favoring certain, more specific priors for the model parameters, 
then the shape of the probability distributions derived from fits will change accordingly, expressing the influence of 
the adopted priors. 
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POLYNOMIAL BASIS FUNCTIONS 



To parameterize supernovae kick distributions, StarTrack employs three parameters, cti, <72, and s, which represent 
the superposition of two MaxwcUian kick distributions with probabilities s and 1— s. The physical predictions associated 
with (cri,(T2.s) arc therefore identical to those of (c72,cri,l — s). To improve the physical significance of our fit, we 
have chosen to employ basis polynomials which enforce this requirement to all orders. 

Specifically, rather than allow for homogeneous basis polynomials in these parameters, we use the following, for 
arbitrary p and q: 



Because ,s must enters in a heterogeneous manner to preserve our desired symmetry, these basis polynomials are of a 
fixed order in all kick parameters. For the purposes of order counting when constructing fits, the first polynomial is 
denoted order p + q and the second order p. 
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